AN INVERSE ITERATION METHOD FOR EIGENVALUE 
PROBLEMS WITH EIGENVECTOR NONLINEARITIES 

ELIAS JARLEBRING*, SIMEN KVAALt, AND WIM MICHIELS* 

Abstract. Consider a symmetric matrix A(v) £ R nxn depending on a vector v £ R n and 
satisfying the property A(av) = A(v) for any a £ R\{0}. We will here study the problem of 
finding (A, v) £ R X R"\{0} such that (A, v) is an eigenpair of the matrix A(v) and we propose a 
generalization of inverse iteration for eigenvalue problems with this type of eigenvector nonlinearity. 
The convergence of the proposed method is studied and several convergence properties are shown 
to be analogous to inverse iteration for standard eigenvalue problems, including local convergence 
properties. The algorithm is also shown to be equivalent to a particular discretization of an associated 
ordinary differential equation, if the shift is chosen in a particular way. The algorithm is adapted 
to a variant of the Schrodinger equation known as the Gross-Pitaevskii equation. We use numerical 
simulations to illustrate the convergence properties, as well as the efficiency of the algorithm and the 
adaption. 



1. Introduction. Let A : R™ — > R" x ™ be a symmetric matrix depending on a 
vector. We will consider the corresponding nonlinear eigenvalue problem where the 
vector- valued parameter of A, denoted v £ R™, equals an eigenvector of the symmetric 
matrix A(v). That is, we wish to find (X,v) elx R™\{0} such that 



(1.1) 



A(v)v = Xv, 



and we will call A an eigenvalue, v an eigenvector and (A, v) an eigenpair of ( 1.1 1. In 
this work we will assume that A is three times continuously differentiable with respect 
to v for any v 0. 

Moreover, we shall assume that A satisfies 



(1.2) 



A(av) = A(v), for any a £ R\{0}, 



such that the solution is independent of the scaling of v, i.e., if (A, v) is a solution 
to then (A, av) is also a solution for any a <E R\{0}. If a problem does not 

have property (1.2), but instead satisfies a normalization constraint, we can often 
transform it to an equivalent problem with property (1.2). For instance, consider 
A(-) which is symmetric with respect to v, i.e., A{v) = A(—v). Also suppose A does 
not satisfy (1.2). The problem to find a normalized vector x € R", i.e., = 1, such 
that A(x)x — Xx is equivalent to (1.1) if we define A(v) :— A(v/\\v\\) such that A 



does satisfy (1.2). 

In this paper we propose a new algorithm for ( |1.1[ ), which is a generalization of 
inverse iteration for standard linear eigenvalue problems. More precisely, we consider 
the iteration, 



(1.3) 



v k +i = a k (J(v k ) - al) 1 v k , 



where a k — l/\\(J(vk) — crI) 1 v k \\2 and J is the Jacobian of the left-hand side of ( 1.1 ) 
with respect to v, i.e., 



(1.4) 



J(v) 



d_ 

dv 



(A(v)v) e 
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Table 1.1 

Applicability of the convergence-characterizations in this paper 



The scalar cr 6 R is called the shift and can be used to control to which eigenvalue 
the iteration converges, and has great influence on the speed of convergence. 



The iteration ( 1.3 ) is a generalization of inverse iteration, which is one of the most 
used algorithms for eigenvalue computations. Convergence analysis can be found in 
[27] H4] and in the studies of Rayleigh quotient iteration in the classical series of 
works of Ostrowski, e.g., [26]. There are also inexact versions of inverse iteration and 
its variant Rayleigh quotient iteration, e.g., [551 IS]- These results are often used in 
combination with preconditioning techniques for inverse iteration [211 [231 E] . Inverse 
iteration has also been generalized to eigenvalue problems with eigenvalue nonlinear- 
ities, e.g., [22l E01 ES] for which convergence has been studied in [15l[T6]. There is 
an algorithm for problems with eigenvector nonlinearities and based on solving lin- 
ear systems in [2T]. See also the results on eigenvector nonlinearities in [7]. To our 
knowledge, the iteration (1.3) for problems of the type (1.1) has not been presented 



or analyzed in the literature. 

We characterize the convergence of the iteration ( 1.3 ) in two ways, which lead to 
conclusions for the behavior in two situations; when the error \\v k — is small or 
when the shift satisfies a <C A*. See Table [TTT] 

In particular, in the local convergence analysis (in Section [3]) we show the follow- 
ing. For any eigenpair (A*,u*) of (1.1), the iteration satisfies 



(1.5) 



v k+1 ± u» = |A» - <r\F*(v k - u») + 0(\\v k - u*|| 2 ), 



where the sign depends on sign(A* — cr), and we derive an expression for 6 R. nxn . 
If the shift is sufficiently close to the eigenvalue, the iteration is locally convergent. 
The convergence is in general linear and the convergence factor is proportional to the 
distance between the shift and the eigenvalue. 



We also provide a characterization of (1.3) by deriving a relation with an ODE 



applicable in the situation when cr < A„. In particular, in Section [4] we show the 
following. One step of (1.3) is equivalent to a particular type of discretization of the 
ODE 

(1-6) y'(t)=p(y(t))y(t)-A(y(t))y(t), 

where p is the Rayleigh quotient 



(1.7) 



p(y) 



V T A{y)y 

y L y 



if the shift is chosen such that a <§C A*. The stationary solutions of (1.6) are solutions 
to (1.1). A small step- length in the discretization corresponds to cr <c A,, implying 



that one step of (1.3) approximates the trajectory of (1.6) if a is chosen sufficiently 
negative. In the linear case (when A(v) = B) when the matrix B has one simple 
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dominant eigenvalue (eigenvalue closest to -co), the ODE only has one stable sta- 
tionary point (which corresponds to the dominant eigenvalue) and it is convergent for 
any starting value. If a similar situation occurs for the nonlinear case, i.e., the ODE 



(1.6) is convergent and only has one stable stationary point, then the iteration (1.3) 
is expected to converge to this solution for sufficiently negative a and the convergence 
is expected to be independent of starting vector. In particular, if the problem is close 
to linear, the iteration is expected to converge to the dominant eigenvalue. 

The idea we use in Section |4j basing the reasoning on interpreting the iterative 
method as a discretization or realization of an ODE, has been used in a number 
of other settings for eigenvalue problems before, e.g., for the QR- method [30] [3TJ 
E], preconditioning techniques [T7j and characterizations of the Rayleigh quotient 



iteration [19]. In fact, the ODE (4.1 ) is a nonlinear variant of the Oja flow [32]. See 
also the collection [3] and the description of iterations on manifolds in [8] . 

The algorithm is illustrated by showing how it can be adapted to the Gross- 
Pitaevskii equation (GPE) in Section [5] The GPE is a standard model for particles 
in the state of matter called the Bose-Einstein condensate. See [TH] and references in 
PQ for literature on the GPE. We discretize the GPE and transform it to the form 



(1.1). We also show how the solution to the linear system (J(vk) — crl) 1 vj c can be 



found efficiently using the Sherman-Morrison- Woodbury formula. It turns out that 



in this setting, the ODE (1.6) is directly related to the technique called imaginary 
time-integration [T]. This connection allows us to use results known for imaginary 
time- integration, in particular that the ODE always converges to a stationary solution 
and the experience presented in results in the literature indicate that this is often a 
physically relevant solution, e.g., the ground state. We further study the ODE and 
derive a heuristic choice of the step-length h, or equivalently a heuristic choice for 
the shift a, which follows the trajectory of the ODE to sufficient accuracy, and still 
maintains a fast asymptotic convergence rate. 



Although we are not aware of algorithms for (1.1), a number of successful algo- 
rithms do exist for the specific application we have in mind, i.e., the GPE in Section[5] 
Besides the aforementioned methods based on imaginary time-integration, there are 
methods for the GPE based on minimization [2j [4] as well as a Newton- Rhapson ap- 
proach [5]. Our approach is different in character since the reasoning stems from an 
eigenvalue algorithm (inverse iteration), and it has a different generality setting. It 
can be interpreted as a discretization of the ODE, but with an integration scheme and 
step-length which we believe has not been used for the imaginary time-integration of 
the GPE. For completeness we also present results for another generalization of the 
algorithm which turns out to be essentially equivalent to an approach presented in [1] 



and can be seen as (1.3) where J is replaced by A. 

The notation in this paper is mostly standard. We use ^q(v) to denote the 
Jacobian of a vector or scalar q, i.e., if v G R. n and q(v) G R k , then ^Lq{v) = 
{j^q(v), . . . , ^-q(v)) G M fcx ". We use the notation { ) x=y to denote substitution of 
x with y for the formula inside the parenthesis. As usual, the expression \\Z\\2 denotes 
the Euclidean norm if Z is a vector and the spectral norm if Z is a matrix. The set of 
eigenvalues of a matrix B G ]R™ X ™ will be denoted \{B) and a dominant eigenvalue will 
be used to refer to an eigenvalue fi G X(B) for which all other eigenvalues have larger 
(or equal) real part, i.e., if \i G X(B) is a dominant eigenvalue then, any /i 2 G A(S) 
satisfies Re (fi) < Re {fii). 

2. Preliminaries and fixed point formulation. Throughout this paper we 
will in several situations use the following consequences of the scaling invariance prop- 
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erty @. 

Lemma 2.1 (Scaling invariance) . Consider a scaling invariant function Q : W 
M. kxe , i.e., Q(av) = Q(v) for any a £ R\{0}. Then, for any vectors v E 



(2.1) 

In particular, given A : I 
defined by (1.4), we have 

(2.2) 



(Q{v)u) )v = 0. 



— >• 



satisfying the scaling invariance (1.2 1 and J 



J(u)u = A{u)i 



for any u E M™. Moreover, ifu — v^, when (A*,f*) is an eigenpair of (1.1 1, then A* 
is an eigenvalue of J(w*) with eigenvector v*. 

Proof. The left-hand side of (2.1) can be interpreted as a directional in the 
direction of v. We have 



(Q(v)u) v 



lim 

e-s-0 



Q(v + ev)u — Q(v)u 



lim 

e->0 



Q((l + e)v)u - Q(v)u 



which shows (2.1). The identity (2.2) follows from the chain rule. □ 
The iteration can naturally be represented in a fixed point form, 

Vk+l = ^(Vfe) 

where 
(2.3) 

and 



\mv)\\ 



(2.4) 



i>{v) := (J(v) - aiy 



Obviously, Vk+i = tp{vk) when is given by (1.3). It also turns out that the fixed 
points of (p are solutions to (1.1). However, the converse is not true. We shall now 
show that if a > A*, we have ip(v*) = — v*, implying that the iterates Vk alternate 
between v* and — v* if vq = v*. An equivalence with the solutions to (1.1) and the 
vectors that satisfy ip(v*) = ±w*, for some choice of the sign, can be achieved if we 
take the alternation into account, as can be seen as follows. 

Proposition 2.2 (Equivalence of fixed points and eigenvectors). 

(a) Suppose (A*,u*) is a solution to (1.1) with \\v*\\ = 1 and suppose a ^ A*. 
Then, the fixed point map ip satisfies 



(2.5) 



The sign in (2.5) should be chosen positive if a < A*, otherwise negative 
(b) Suppose = ip(v*) or u* ■ 
with 



(2.6) 



A* := <7 ± 



</?(?;*). Then, (A*,i>*) is a solution to (1.1) 
1 



\\{J{v*)-o-I)-H4 



A(v*)v* 



•i 



Proof. In order to show (a), suppose (A*,u*) is an eigenpair and note that from 
Lemma [2~T1 we have that 

X^v* = A(v*)v* — J(v«)u». 

By subtracting o~v* from both sides and subsequently multiplying both sides from the 
left by A \ (J(v*) — crl)^ 1 we can simplify 



and we consequently find that 
(2.7) ' 



1 



Hence, </? evaluated at v = is explicitly 

1 



= I A. 



y/^(v t ) T i)(v t ) 



\X* - tr| 

A* — (T 



±u*. 



In order to show (b), suppose that is such that </j(i>*) = ±t!«. Equation (2.3) leads 
to 

(J(u*) - o-I)u* = n _ ln ^- 
||(J(i;*) — o-J) !|| 



From Lemma 2.1 it follows that J(v*)v* — A(v<,)vt, and we finally have that 

1 



A(v*)vi, = [a ± 



IKJ^-a/)- 1 ^!! 



and the formula (2.6) follows by multiplying from the left with v£ . □ 
3. Local convergence properties. 

3.1. First-order behavior and convergence factor. From Proposition |2.2| 
we can directly conclude that if we at some point in the iteration have Vk = i>*, where 
v,,, is an eigenvector, then every subsequent iterate Vj, j > k, will also be an eigenvector 
corresponding to the same eigenvalue, but can possibly alternate between ±u*. We 
will now study the case where Vf. is close but not equal to an eigenvector. In order to 
understand this local convergence behavior, we also need to take the alternation into 
account. The error is to first order given by the following result. 

Theorem 3.1 (Local convergence). Suppose (A*,u*) is a solution to (1.1) with 
Hu* || = 1. Let a € K be any shift such that J(v*) — a I is non-singular, in particular 
a 9^ A* . Then, 



(3.1) 



<p'(y*) — | A* — a\(I — v*vj) (J(i>*) — o-iy 



Moreover, suppose the iterates Vk generated by (1.3) are such that J(vk) — al is non- 
singular for any k. Then, the iterates satisfy 



(3.2) 



Ufe+i T v* = ^>'{v*)(vk - «*) + 0(\\v k - f*|| 2 ) 



The sign in (3.2) should be chosen negative if a < A*, otherwise positive. 
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Proof. By using the Taylor expansion and the fixed point characterization in 
Proposition 2.2 we have that the left-hand side of (3.2 1 satisfies 

v k+ i =f v* = v k +i - (±««) = <p(vk) - ¥>(«*) = - "*) + - w*|| 2 ). 



It remains to derive the formula ^'(f*) given by (3.1 1. 

We will need the derivative of the two-norm, which can be expressed as the row- 
vector 



(3.3) ^\mv)h = ^m T m = 



1 



y/ip{v) T ^(v) 



^{vYi)'{v) = ip(v) l ip'{v) 



since f{v) T — ip(v) T /\\ip(v)\\2- 

From the definition (2.3) of ip we have the relation y>(?;)||'0(t;)|| = i/j(v), whose 



Jacobian can be computed with the product rule and ( |3.3[ ), 

<p'(v)\\i>(v)\\ + tp(v)tp(v) T ip'(v) = ip'(v) 

Hence, 

1 



(3.4) 



P» = 



\Hv)\\ 



(j-vj(«M«) t )^(«). 



Now recall (from Proposition 2.2) that <p(v*) = ±v* and consequently ip(y tf )tf(v ie ) T = 
Vstvl '. This fact combined with the formula for the norm (2.7) leads to a simplification 
of (|3.4|) when we evaluate at v = u*, 



(3.5) <f/(v m ) = | A, -a\(l- v m vl) t/>>*). 

It remains to establish a formula for the Jacobian of ip evaluated at v — u*. By 
differentiation of (2.4) multiplied by J(v) — al, we have 

(3.6) 



-(J(v) - al)ip(v) ) + (J(v) - (rl)ij'(v) 

V I - 

/ v— v 



and 



(3.7) ip'(v*) = -(J(v*)-aiy 



J(v)ip(v* 



+ (J(v t )-aiy 



Moreover, note that ip{v*) — || ip(v* )|| tp(v*) — ±||(J(u*) — al)~ 1 v*\\v*. We will now 
show that the first term in ( |3.7[ ) vanishes identically by showing that all columns of 
(,~§ijJ( v ) v *)v=v vamsn - Let the jth column be denoted 



(3.8) 



J(v)v t 



ej = lim - (J(u* + eej)v sr — J(y*)v M ) . 



Lemma |2.l| and in particular the identity J(u)u = A(u)u for any u £ R™, implies that 

J(w* + eej)v* — J(v»)u* = 

J(w* + £ej)(u* + ecj) — eJ(v* + sej)ej — J(f*)u* = 
+ sej)(v« + eej) — eJ(v« + sej)ej — J(v*)v* = 

j4(v*)u* + eJ(v*)ej — eJ(u* + sej)ej — J(v*)u* + 0(e 2 ), 
6 



where we used a Taylor expansion of A(v* + eej)(v* + ee 7 ) = A (v*)v* + J(v*)(eej) + 
0(e 2 ) in the last step. Hence, by again applying Lemma 2.1 now giving A(v*)v* = 
J(u*)v*, we have 

cj = lim (j(u*)ej — J(u* + £ej)ej + 0(e)) = 0. 
We have shown that Cj = 0, for any j = 1, . . . , n which implies that 

(3.9) J{v)v* J = e M" x " 



such that the first term in (3.7) vanishes and the proof is complete. □ 

The above theorem for the local behavior directly gives a characterization of 
the convergence factor in terms of the largest eigenvalue of tp'{v*) in modulus. To 
formalize this statement we will use the concept of Q~convergence factor and R- 
convergence factor given in |25[ Chapter 9]. We briefly summarize the concepts for 
the generic situation and our setting. Suppose a sequence {vfcj^Q converges linearly 
to u*. Then, generically, the R- convergence factor (for the sequence {vk}%^ =0 ) is given 

by 



Ri({v k }) = lim \\v k - 

k— foo 

Correspondingly, the Q-convergence factor is given by 

||ufe+i - v* 



ii/fc 



Qi ({«*})= lim 

fc->oo \\v k — V*\\ 

The i?-convergence factor and Q-convergence factor associated with a fixed point it- 
eration and a particular fixed point is defined as the supremum of the convergence 
factor for all non-degenerate convergent sequences generated by the fixed point iter- 
ation, and will be denoted R\ (</5,f*) and Qi((p,v*). See [2SJ Chapter 9] for formal 
definitions. 

When a < A* we can conclude from Proposition |2.2| that is a fixed point of 
<p(v*). On the other hand, if a > A* we need to take alternation into account. It is 
easy to verify that is fixed point of <p{v) := —<f(v) which is a fixed point iteration 
giving the same sequence of vectors v k generated by ip, except that every second vector 
has a different sign. This compensates for the asymptotic sign alternation, such that 
convergence to an oscillating sequence v k using tp is equivalent to a convergence to a 
fixed point using (p and allows us to study both cases in a unified manner. We will now 
see that these definitions of convergence factors give one formula for the convergence 
factor (independent of the two cases A* < a or A* > a) involving eigenvalues of J(w*). 
Moreover, the i?-convergence factors and Q-convergence factors are in the generic 
situation equal. 



Corollary 3.2 (Convergence factor). Let (A*,u*) be an eigenpair of (1.1) with 
||i;*|| = 1. Suppose A* is a simple eigenvalue of J(y*) and let 

(3.10) 7 -l A *- ff l 



|M2 - q 

where //2 S C is the eigenvalue of J(««) closest to a, but not equal to A*, i.e., 

/i 2 = argmin | — cr | 

ve\(J(v,))\{\,} 



and suppose it is simple. Let <p be defined by 
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• ip := ip when a < A*; and 

• ip := — u;/ien er > A*. 

Suppose 7 < 1. Then the fixed point v* is locally convergent with respect to ip and the 
convergence is linear (in both R order and Q order). Moreover, the R- convergence 
factor is given by 

(3.11) fl 1 (^,«0=7- 

Furthermore, assume that {wfe}^ is a sequence generated by ip convergent to 
and such that R\({vk}) = 7- Then, under the condition that \\vk+i — — 
converges to a nonzero value, the Q-convergence factor is equal to the R-convergence 
factor, i.e., 

(3-12) Qi({%})=7- 



Proof. Consider any of the two iterations Vk+i = <p(vk) and ffc+i = —ip(vk). 
Ostrowski's theorem [25J Theorem 10.1.3], the linear convergence theorem [25], The- 



orem 10.1.4] and Theorem 3.1 gives the characterization of the i?-convergence factor 
using the spectral radius 

fli ({«*}) = P(V - v*v?)(J(v.) - al)- 1 ), 

independent of which of the two iterations are considered. 

Now note that is an eigenvector of J(v*) and also an eigenvector of (J(v*) — 
cr/) _1 . The application of the projector (I — v^vj) has the effect that the eigenvalues 
of (I— v*v T )(J(v*) — o-I)^ 1 are the same as the eigenvalues of (J(v^)-al)^ 1 except for 
the eigenvalue corresponding to (which by assumption is simple). This eigenvalue 
is transformed to zero, hence the maximum must be taken over the eigenvalues of 
(J(u*) — o-I)^ 1 except for the (er — A*) -1 . We shown ( 3.10[ ) have the formula for 7. 



The formula (3.12| follows from the fact that the linear Q-convergence factor is 
equal to the linear i?-convergence factor under the assumption that the error quotient 
||Vfc+i — w*||/||wfe — converges, which holds by assumption. See e.g. [231 NR 10.1-5]. 
□ 

3.2. Similarities and differences with inverse iteration for linear prob- 



lems. The iteration (1.3 1 is clearly a generalization of inverse iteration when A(v) — 
A is constant, and a behavior similar to (linear) inverse iteration is expected when 
the problem is close to linear. More importantly, from the theory above we can con- 
clude that the iteration possesses many properties similar to inverse iteration also in 
a situation when the problem is not close to linear. 

• The convergence is in general linear. 

• The convergence factor is asymptotically proportional to eigenvalue shift dis- 
tance when the distance is small. 



The convergence factor (3.10) is in general determined by a quotient involving 
the eigenvalue shift distance in the numerator. Unlike the linear case, the de- 
nominator is the distance between the shift and the second closest eigenvalue 
of J(w*). 

Unlike the linear case, for a given shift, several eigenvectors can be attractive 
fixed points. Conversely, for a given shift, the iteration may not have any 
attraction points. 



In order to further illustrate the value of these properties we will in this work also 
briefly study another generalization of inverse iteration. Inverse iteration for standard 
eigenvalue problems consists of shifting, inverting and normalizing, and it is from this 
perspective natural to consider the generalization Vk+i = <PA(vk), 



(3.13) 



<Pa{v) ■-- 



We will call the A-version of inverse iteration. 

We can derive a description of the first-order behavior similar to Theorem |3.1| 
Note that ip' A (v*) is considerably more complicated than <f/(v^] and most of the bullets 
above do not apply to this version, most importantly, the convergence factor is not 
necessarily small if the shift is close to the eigenvalue. 

3.3. Illustration of local convergence. Several properties of the algorithm, 
including the local convergence above, can be observed when applied to the example 



Bv 

A[v) = A + sin( — = — )^!. 



The Jacobian is given by 



J(v) = i-A(v)v = A(v) + 2 . ' A lV ((v T v)v T B - {v T Bv)v T ). 
av (v 1 v) z 

We selected Aq, A\ and B in a random way. In order to make the numerical simula- 
tions reproducible, we will fix Aq, A\ and B as follows 
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where we carry out simulations for a number of different /3. Note that the problem is 
linear when {3 = 0. 

Simulations of the inverse iterations for different /3 are given in Figure |3.1| and 



Figure 3.2 In Figure 3.1 we clearly observe linear convergence, for the A- version 
as well as the J-version. The J-version is also clearly less dependent on /?. As the 
nonlinearity parameter f3 increases, the performance of the A-version worsens whereas 
the convergence of the J-version is not greatly affected by increase of j3. We see in 
Figure [3~Tf that the slow convergence of the A- version can not always be compensated 
by moving a closer to the eigenvalue A*. In this case the A- version is not even locally 
convergent when a — A*, whereas the J-version exhibits quadratic convergence. Note 
that for a = A* the fixed points map ip(v*) and lpa{v*) involve inverses of singular 
matrices, implying that they do formally not have fixed points v*. is invertible. 
Similar to inverse iteration for standard eigenvalue problems the J-version still works 
in practice when subject to rounding errors. 
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(a) p = 0.01, cr = A* +0.3 



(b) /3 = 0.1, <r = A, + 0.3 
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(c) P = 0.5, cr = A, + 0.3 



(d) P = 0.6, cr = A, + 0.3 
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(e) p = 1, cr = A» +0.3 



(f) P = 1, a = X, 



Figure 3.1. Convergence for the example in Section \3.3\ using \1.3\ , i.e., J-version of inverse 
iteration, and the A-version of inverse iteration l|3.13|l. 



In Figure 3.2 we see that the convergence factor of the J-version approaches zero 



when a — > A* as predicted by Corollary 3.2 The convergence factor for the A-version 
does clearly not vanish when a is close to A* . 



4. Interpretation as the discretization of an ODE. 
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(a) 



(b) 



Figure 3.2. Estimated convergence factor for J -version and A-version of iteration, where 
estimation is done with the quotient ||iifc+i i vk ||/||"fc ± "if || where vk is a very accurate solution. 
Clearly, when a — > A* the convergence factor for the J -version approaches zero, unlike the A-version. 
When a — > — oo the convergence factors coincide. (5 = 0.5. 



4.1. The normalized ODE associated with In order to provide fur- 

ther insight into the non-local behavior of the iteration, we will consider the following 
ODE 

(4.1) z'(t) = -A(z(t))z(t) 

and the function corresponding to the normalization of the solution 

(4.2) y(t) := q(t)z(t) where q(t) 



1 



\z(t)W 

The function y also satisfies an ODE, which does not involve z. This can be seen from 
the following reasoning. Note that 



2z(t) T z'{t) = 



2(z(t) T z(t)) 3 /< 



1 



:(t) T A(z(t))z(t) = q (t)y(tfA(y(t))y(t). 



By applying the product rule to (4.2), we have 

y'(t) = q(t)z'(t) + q'(t)z(t) - - q(t)A(z(t))z(t) + q(t)(y(t) T A(y(t))y(t))z(t). 
Hence, the normalized function y(t) satisfies the differential equation 
(4-3) y'(t)=p(y(t))y(t)-A(y(t))y(t), 



where p(y(t)) is the Rayleigh quotient defined by (1.7 1. The ODE ( |4.3[ ) will be used 
as a tool to understand the iteration (1.3). For this reason several properties of the 
ODE and some of its relations with the nonlinear eigenvalue problem ( |1.1[ ) will be 
particularly important. 



Theorem 4.1 (ODE properties). Consider the ODE (4.3) with initial condition 
y(0) = yo with \\ya\\ = 1. The ODE and its solution have the following properties. 
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(a) The norm \\y{t)\\ = 1 independent oft. 



(b) Any stationary point of the ODE (4.3) is a normalized eigenvector of (1.1) 



(c) Any normalized eigenvector of (1.1) is a stationary point of the ODE (4.3) 



(d) Let y* be a stationary point of the ODE (4.3). The stationary point is stable 
if the eigenvalue A* = p{y*) is a simple dominant eigenvalue of J{y*), i.e., if 
A* is simple eigenvalue of J{y*) and 

A* < Re (/i) for any fi € A(J(y*))\{A*}. 



(e) Let y* be a stationary point of the ODE ( 4.3 I . The stationary point is unstable 



if there is an eigenvalue (i of J(y*) such that 

Re (p) < A* =p(y„). 



Proof. The statements (a)-(c) follow directly from the derivation of (4.3). In 



order to study the stability of the stationary point g/*, let A(t) := y(t) — y* denote the 



deviation from the stationary solution. Then, by Taylor expansion of (4.3) we have 

A'(t) = y'{t) = p(y(t))y(t)-A(y(t))y(t) = \p(y*)I + y* P '(y*) - J(y*)} A(t)+0(A(t)) 2 . 

We can now insert p(y*) = A* and p'(y*) — yJ{J(y*) — A/) since = \\yo\\ = 1. 
Hence, the behavior is to first order given by 



(4.4) 
where 



A'(t) = F(y t )A(t), 



F(y*) = X,I + y*yJ(J(2/*) - XI) - J(y m ) = (/- y*yl)(\J - J{y*)) 



The eigenvalues of F(y*) have the form 



A* 



where /x is an eigenvalue of J{y*), but not equal to A*. The ODE (4.4) is unstable 



if F(y*) has eigenvalues with positive real part, and the ODE (4.4) is stable if all 



eigenvalues have negative real part. The statements (d)-(e) follow from the sign of 
Re (A* — fi) = A* — Re (/i), which, provided it is not zero, gives a conclusion about 
the local stability of (O) and p~3| ). 
□ 

4.2. Discretization of the ODE and equivalence with inverse iteration. 



In the previous subsection we saw how the ODE (4.3) was related to the nonlinear 



eigenvalue problem (1.1); in particular we showed that the stationary solutions were 



equivalent to the eigenvectors of (1.1). We will now show how this connection can be 



used by showing that the proposed version of inverse iteration ( 1.3 ) can be interpreted 



as a discretization of the ODE, allowing us to study the general behavior of the 
iteration by studying the ODE. 

Consider first the Rosenbrock-Euler method O Chapter IV. 7], i.e., the backward 
Euler method where the linear system is solved with one step of Newton-Raphsons 



method. The backward Euler method with step-length h applied to (4.3) is 



(4.5) 



Vk+i ~ Vk 



p{yk+i)yk+i 
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A{vk+\)yk+i- 



We approximate the first term in (4.5) by its linearization, 



p{yk+i)yk+i ~ p{yk)yk + (p'(yk)yk + p(yk))(yk+i - Vk) = p(yk)yk+i, 

where we denoted p'(z) := -§zp{z) & M lx ™. We also used that p'(z)z = for any 



z which follows from Lemma 2.1 since p(z) is scaling invariant. The second term in 
(4.5) is now also approximated by its linearization. By combining the approximation 



with ( 2.2 ), we have 



Mijk)yk + J(vk)(Vk+i - Vk) = J{yk)yk+i- 

Hence, the Rosenbrock-Euler method applied to ( |4.3| ) is equivalent to 

frWk+i - Vk) = p{yk)yk+\ - J{yk)yk+i 

and also 

(4-6) ((l-p(y k ))l + J(y k ) 



1 



Vk+i = vVk- 



The ODE (4.3) has, according to Theorem 4.1 a solution with constant norm. 
Despite this, the discretization y k does not have constant norm, due to the approx- 
imation error introduced in the stepping scheme. This issue can be addressed by 
carrying out a standard projection described, e.g., in [lit Algorithm IV. 4. 2], which is 
a common procedure to incorporate normalization constraints. In this case, it reduces 
to subsequently normalizing the iterate, i.e., setting y k +\ '■= jrr^ — nVk+i- Since p and 

II Vk + l || 

J are scaling invariant, we have p(jjk) — p{yk) and J(yk) = J(yk), and we can express 
one step of the integration scheme as 



(4.7) 



Vk + l = OL k 



h ~ p ( yk "> 



J{Vk) 



Vk, 



where a k = 1/ ((£ - p(yk)) I + J{Vk))~ 



Vk 



The main result in this section is based on the observation that (4.7) is equivalent 



to the proposed version of inverse iteration (1.3) if the shift a and step-length h are 



related in a particular way. This connection is summarized in the following theorem. 

Theorem 4.2 (ODE discretization equivalence). Let the sequences {j/fc}^L and 
{"^fclfeLo ^ e generated in the following way. 



(a) Let v k be the iterates generated by inverse iteration with (1.3) with a given 
shift a e E. 



(b) Let y k be the discretization of (4.3) using 



the Rosenbrock-Euler method, 

a standard projection step which imposes the normalization, and 
the step-length 

(4.8) I ' 



p(Vk) - o- 

where a is the shift used to generate v k . 
Moreover, suppose they are started in the same way, i.e. 

v k = y k 

for all k £ N. 

We summarize some immediate consequences. 
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Figure 4.1. Convergence of the ODE (sine- example) to the dominant eigenvector. 



The iteration (1.3) is an approximation of the trajectory of the ODE (4.3) if 
h is chosen sufficiently small, i.e., if a is chosen sufficiently negative. 
Suppose the ODE converges to a stationary solution, which indeed is the 
case for several examples, and can in particular be proven for the example in 
Section [5j In this situation the stationary solution is a solution to the ODE, 
and the iteration will converge to an eigenvector for sufficiently negative a . 
For the linear case, i.e., A(v) = B, we can explicitly write down the solution 
to (4.1 ) using the Jordan form, and thereby study (4.3) which will converge to 



the dominant eigenvalue unless y(0) is chosen in a particular way. Moreover, 
Theorem |4.1| i-e shows that if the dominant eigenvalue is simple, it will be 
the only stable solution. Analogously, for nonlinear problems which resemble 



the linear case in the sense that the ODE (4.3) only has one stable stationary 
solution and the ODE is convergent (unless started in a particular way), the 
iteration will converge to that solution for sufficiently negative a. 
Remark 4.3 (ODE interpration for A-version). The A-version (3.13) can be 



interpreted as discretization of the ODE (4.3), but where J is approximated by A 



and h is chosen according to (4.8). In fact, when the A-version is applied to the 



example in Section [5] it is equivalent to an available in the literature; more precisely 
the discretization of the normalized gradient flow given in JIJ Equation (2.20)-(2.21)]. 

4.3. Illustration and numerical justification of ODE discretization. Con- 
sider a situation where we do not have any approximation of any eigenvalue of the 



example in Section 3.3 To characterize this situation we carry out a number of simu- 
lations with random starting vectors for different shifts. Such a statistical simulation 
is shown in Figure |4.1fc ,. We observe that the for a < —5, the iteration appears to 
always converge to the dominant solution. This property can be explained using the 
developed theory. When a <C A* we can apply Theorem 4.2 and when \o~ — A*| <C 1 
we can apply the local convergence theorem Theorem 3.1 



We computed the solution for the ODE (4.3) for this example using a standard 
ODE integrator, with a random starting value. The trajectory is shown in Figure [4~Tp - 
c. Clearly, the ODE converges to the dominant eigenvector. Similarly, in Figure |4.2| 



we have shown the computed approximate trajectories using (1.31 (J- version) and 
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(3.13) (A- version). We see in accordance with the interpretation of a as a particular 



choice of step-length (via ( 4.8 )) that the iteration follows the trajectories better if a 
is more negative. Moreover, with the J-version we can take larger steps than with 



the A-version, as the A-version does not follow the ODE in Figure 4.2 1 



The experience with the ODE (1.3) for this particular example indicates that the 



ODE converges to a stationary solution for all /3 e [0, 1]. Moreover, for several choices 
of /?, among all eigenvectors, there appears to be only one situation where A* is the 
dominant eigenvalue of J(u*). This can be seen in Figure 4.3 From Theorem 4.1 



we consequently know that the ODE (1.3) only has one stable stationary solution 



Combined with the convergence observation, the ODE will always converge to this 
solution. Hence, Theorem |4.2| implies that for sufficiently negative a, we will follow 
the trajectory of the ODE and eventually converge to the dominant eigenvector. 




(a) cr = -7 





0.5 1 1.5 2 




(b) a = -10 



(c) a = -20 



(d) a = -50 



Figure 4.2. J-version is better in the sense that it follows the ODE also for larger h than for 
the A-version. The dominant eigenvalue is A* Ri —6.01. 



5. Application to the Gross-Pitaevskii equation. 

5.1. Discretization and reformulation to the form To a achieve 

the physical phenomeon of Bose-Einstein condensation, weakly interacting particles 
(such as 4 Hc atoms) are trapped and cooled down to very low temperatures, thereby 
undergoing a phase transition into a Bose-Einstein condensate (BEC), a superfluid 
phase of collective quantum mechanical motion. A standard model for a BEC is the 
Gross-Pitaevskii equation (CPE), a nonlinear partial differential eigenvalue equation. 
This equation exhibits interesting topological behavior. In particular, it reproduces 
experimentally observed quantized vortices in a BEC within a rotating trap, the 
hallmark of superfluidity. 
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Figure 4.3. The eigenvalues of J(t>*) for the example discussed in Section \4-3\ In every 
subfigure, □ denotes the eigenvalues of J(v„) for an eigenpair (A*,ti„) of and A, is denoted 

by X. In this notation, a stable stationary solution described by Theorem \4- l\ l-e corresponds to a 
situation where K is to the left of all □. 



Upon discretization, the GPE becomes an eigenvalue problem of the form (1.1) 



and to illustrate the iteration \l.2>\ we will consider the case of a rotating BEC on the 

domain D = (— L, L) x (— L, L) in two dimensions, whose GPE reads 

(5.1) 

Id \ 

■-A-ifi— + V(x,y) ip{x,y) + b\tp(x,y)\ 2 ip(x,y) = Xi(j(x,y), (x,y) G D, 



where ip(x, y) is the condensate wavefunction (not to be confused with ip in Eqn. (2.4)) 



such that \ip(x,y)\ 2 is proportional to the particle density at (x,y), and where 5 is 
the angular derivative given by 



(5.2) 



^ dx dy 



The function V(x, y) is the external potential function which describes the shape 
of the particle trap. Here, we choose an asymmetric harmonic oscillator potential 
V(x,y) — (x 2 + \.2y 2 )/2. We are mostly interested in those solutions to (5.1) which 
are physically important, e.g., the ground state. 

The boundary condition is chosen as ip{x,y) = on (x, y) G <9B. Moreover, tj) 
is normalized such that ||V'IIl 2 (d) = 1- The physical constant b controls the strength 
of the interactions between the bosons, and £1 the angular velocity of the rotation, 
which are here selected to b — 200 and £1 = 0.85, respectively. 

The domain is discretized with N + 2 equidistant grid points in each physical 
direction leading to n c = iV 2 interior grid points with spacing Ax in both directions. 
For completeness we provide the matrices resulting from the discretization. We use 
the approximations of A and d/d(f>, 



Ln = D 2 ,n®I + I®D 2i n, L^.n = diag(yi, . . . , y N ) ® D N - D N ® diag(^i, . . . ,x N ), 
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where Dn and D 2 ,n are the central difference approximation of the derivative and 
the second derivative respectively. Moreover, we let 



V = (V(xx,yx), . . . ,V(x N ,y 1 ),V(x 1 ,y 2 ), . ■ ■ ,V(x N ,y 2 ), . . .V(x N ,y N )) e R n 
Then, with 

A = ~2 L N - ittL<p >N + diag(V r ) 

the discretized problem is 

(5.3) A{z)z = Xz 
where 

(5.4) A(z) -i +/3diag(|z|) 2 , 

and to be consistent with ||"0||l 2 (o) = 1 we must have %p(xj,yk) = (Aa;) -1 zjv(fc-i)+j 
and /3 = (Ax)~ 2 b. 

The problem ( 5.3 1 is a complex problem of dimension n c , not satisfying the scaling 
invariance (1.2 1. In order to transform this problem to the form ( |1.1[ ) we introduce 



the following notation 
A(v) :-- 



ReA( - ) -Im A( - v ^ +m \ 

Im A( . p±S= ) Rei( , Z 1+ ™ 2 T ) 



Re i -Im A \ P , , 
Imio Re i / « T w l J ' 

and 

n,. ._ Mag(wi) 2 +diag(w 2 ) 2 \ _ /S 

W "~ ^ diag^) 2 +diag(t> 2 )V ' V ~ [v 2 

With this definition of the matrix A(v) we have transformed A{z) by treating the real 
and imaginary parts of z as separate variables. Consequently, the complex problem 



(5.3) of dimension n c is equivalent to the real problem 



(5.5) A(v)v = Xv, 



of dimension n — 2n c which does satisfy (1.2) 



5.2. Jacobian and exploitation of Jacobian structure. The resulting ma- 



trix A(v) in (5.5 1 is a sparse matrix for any v due to the finite-difference discretization. 



In order to carry out (1.3) we need the Jacobian J(v) and we need to be able to ef- 
ficiently solve the corresponding shifted linear system of equations. It turns out that 
the Jacobian is in general a full matrix, making the direct application of sparse solvers 
inefficient. Fortunately, the linear system involving the Jacobian can be decomposed 
into two linear systems of equations involving a matrix which is sparse. The deriva- 
tion is based on the fact that the Jacobian is the sum of a sparse matrix and a 
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rank-one matrix. For such structures the Sherman-Morrison- Woodbury formula [101 
Section 2.1.3] is a standard technique. 

Theorem 5.1 (The Jacobian associated with the Gross-Pitaevskii equation). 
The Jacobian for the nonlinear eigenvalue problem (5.5 1 corresponding to the Gross- 
Pitaevskii equation is given by 



(5.6) J(v) 



8 



do 



(A(v)v) 



/Re A Q -Im A 

\lm A Re A 

3diag(«i) 2 + diag(v 2 ) 2 2 diag(ui) diag(v 2 ) 

2diag(«i) diag(v 2 ) diag(^i) 2 + 3diag(u 2 ) 2 



where v T = (vj r ,v 2 r ) with V\,1)2 & M™/ 2 . Moreover, suppose v T v 



- -7f^B(v)vv T 
v 1 v 

1 and let 



C := 

Re A -Im A 
Im A Re A 



3 diag(wi) 2 + diag(w 2 ) 2 2 diag(wi) diag(w 2 ) 
2diag(wi) diag(w 2 ) diag(^i) 2 + 3diag(v 2 ) 2 



- al 



and suppose C is non-singular. Then, the solution to the linear system (J(v)-crl) 1 v 
can be expressed as 



(5.7) 

where 



(J(v)-crl) v = u\+U2, 
ui = C~ x v 



U2 



T 
V Ui 



w with w = 2/3C~ 1 B{v)v. 



Proof. The formula for the Jacobian (5.6 1 follows from several applications of the 
product rule. More precisely, we have 



8 



8v 



A(v)v 



8 



v- 1 v 
(v T v) 2 



Y^B{v)v 



T 

V V 



'/Re A) -Im A 
dv V V Im A o Re 
'Re A -Im A " 
. Im A Re A 

The term involving the Jacobian of B(v)v can now be simplified, 
d 



d_ 

dv 



B(v)v - 2B(v)w 1 



,-B(v)v 
ov 

3diag(u) 2 + 



0_ 
dv 

I 

1 



diag(u) 2 « + 
diag(w) 2 



I 

1 



I 

1 



diag(v) 2 
2 



B(t;) + 2diag(t;) 2 + 2 



I 

1 

I 

1 



I 

1 



diag(w) 
diag(w) 



I 

1 

I 

1 



diag(u) 



diag(w) 



I 

1 

I 

1 



The relation (5.6) follows from expanding and combining the three terms. 

The formula for (15. 71 follows from the fact that the last term in (15. 61) is a rank-one 



term and we can apply the Sherman-Morrison- Woodbury formula |10i Section 2.1.3]. 
□ 
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5.3. Specialized ODE interpretation and step-length heuristics. We will 



now use the following important observation. The ODE (4.3) is equivalent to the 



ODE representing a flow in the technique known as imaginary time propagation or 



normalized gradient flow in, e.g. [Tj. In particular, (4.3) is equivalent to [H Equation 



(2.15)-(2.16)]. This connection allows us to directly reach conclusions about the ODE 



(4.3) for the GPE. We conclude from (TJ Remark 2.6] that the ODE will converge 
to a stationary solution. It can equivalently be shown to converge by constructing a 
Lyapunov function and applying a variant of Lyapunov's second method. (See [131 
Chapter 3] for Lyapunov's second method.) Moreover, the imaginary time propa- 
gation technique is for physical reasons known to converge to a physically relevant 
solution (there may be several), e.g., the ground state or meta-stable configurations 
of the system. 

In light of the equivalence with imaginary time propagation, it is natural to follow 
the true flow as closely as possible during the iteration, while at the same time taking 
as long time steps as possible to minimize computation time. Therefore, we propose 
the following heuristic for the step length h or equivalently a choice of a. 

Suppose that we have computed an approximation y k m y(t k ), and we now need 
to determine an appropriate h to compute y k +\ ~ y{tk + h). A step-length choice 
for the Rosenbrock-Euler method is given in Appendix [X] and in particular formula 



(A. 2), for a given fixed local error e. Note that for our case, i.e., the ODE (4.3) 



f(v) = p(v)v — A(v)v and 

f(v) = -{I - vv T )(J(v) - p(v)I) + vv T {A(v) - p(v)I) 

such that 

f'(v)f(v) = (I- vv T )(-J(v)f{v) + p(v)f(v)) + vv T (A(v) - p(v)I)f(v). 

We will now see that this expression can be efficiently computed before each step of 
the iterative method. 

This leads us to the choice of the shift a we propose to use in this work: 



Compute fk 
Compute <7fc 
Compute e k 
Compute hk 



= /(«*) = p{vk)Vk ~ A(v k )y k 
= J(vk)fk 

= (I - v k vl)(-e k +p(v k )f k ) + v k vl(A(v k ) -p(vk)I)fk 
= (2e/\\e k \\) 1 / 2 using a desired tolerance e. 



If hk > /imax set h k = h 



max • 



• Set a = p(v k ) — 1/hk according to (4.8). 
The choice to make sure that h k < h max is to avoid taking too large steps, for which 
the reasoning for step-length above is not supported. 

Note that we do not need form the matrices / — v k v^ or J(v k ) explicitly when we 
apply it to the GPE. It is more efficient to instead compute the vectors f k , g k and e k 
by forming products between vectors and matrix vector multiplications with sparse 
matrices, since, e.g., J(v k ) is given as the sum of a matrix and a rank one matrix in 



Theorem 5.1 1 



5.4. Conclusions from computational results. We carried out the inverse it- 



eration algorithm with the heuristic choice of a proposed in Section 5.3 for a number 



of different choices of parameters. Wc selected b = 200, L — 15, and 51 = 0.85. 



The number of grid points is N = 300, i.e., the eigenvalue problem (1.1) is of 
size n = 180000. The step-length heuristic was chosen with parameter e = 2 and 
^max = 10 4 (except where otherwise stated). An initial guess was chosen as a random 
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superposition of Gaussians, such that its length scale is independent of the interior 
grid size N. The simulation was completed in 1.2 • 10 3 seconds with an implemen- 
tation of the algorithm in MATLAB running on an Apple MacBook Pro with a 2.6 
GHz Intel i7 quad-core processor. 

The results of the numerical simulations are presented in Figure 5.1|5.5 In these 
figures ipk denotes the approximate solution to (5.1) after iteration k iterations, and 
K denotes the total number of iterations. 



The convergence is visualized in Figure [572] As expected from the ODE interpre- 
tation in Section [4] and the fact that the GPE ODE converges, we eventually reach 
a solution. Moreover, the asymptotic convergence is fast as the step-length is larger 
when the solution is accurate. The approximations of the solution at different iterates 
are given in Figure 5.3 showing the shape of the function as it evolves and converges. 
Clearly, the random initial condition turns into a physically meaningful approxima- 
tion after only a few iterations, and the final iterations mostly change the position of 
the vortices. 

We can also observe the local convergence properties presented in Section [3] in 
this application. Note first that when h — /i max , we can see that a is almost constant. 
This can be observed in Figure |5.5| It is also expected from the fact that the error 
behaves like \\vk — vk\\ ~ \\y{tk) — y*\\ oc C exp(— atk) in an asymptotic regime. 
This can indeed be confirmed in Figure |5.2| The estimate clearly indicates that 
\\ v k+i ~ ^A'll/ll^fe — v k\\ should converge when tk+% = tf. + /i max . Hence, with the 
assumption that a is approximately constant in the regime where we take step-length 
/[ max we have a rj A* — l/h max . We can compute the theoretical convergence factor 
for this a by using Corollary |3.2| and computing [ii with the Arnoldi method (and 
a matrix- vector product from Theorem 5.1|. The theoretical convergence factor is 
visualized together with the estimated convergence factor ||f/c+i — ^i<-||/||i>fc — vk\\ in 
Figure |5.4| The theoretical convergence factor is confirmed for two different choices 

of h ill;! \' • 



The heuristic choice of a is visualized in Figure 5.4 With the crude estimation 



of the relation (4.8), h ~ 1/(A* — c), we see that the step-length in the beginning is 
chosen large, in an intermediate phase it is chosen around the order of magnitude 10 
and in the final phase it is chosen larger, and a is again eventually almost constant. 



6. Concluding remarks. The favorable convergence properties of the cele- 
brated inverse iteration algorithm for the standard eigenvalue problem are well under- 
stood. An important point in this paper is that the generalization we have presented 
have many of the favorable properties that are present in the inverse iteration algo- 
rithm for standard eigenvalue problems. This holds in particular for local convergence 
and the interpretation as an ODE. We have also illustrated the usefulness of the al- 
gorithm by adapting it to a variant of the Schrodinger equation. 

The connection between the GPE and the use of inverse iteration presented in 
this paper has further indirect value. For instance, the tremendous amount of under- 
standing that is available for inverse iteration (for standard eigenvalue problem) now 
have the potential to be exploited or adapted to this type of nonlincarity. 
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Figure 5.1. Left: Visualization of the paricle density \ipK{x, y)\ 2 of the computed solution. 
The contour levels are selected as 0.05 ■ 10 -3 , 0.10 ■ 10 -3 , 0.30 ■ 10 -3 , the outer contours 

have smaller values. Right: Visualization colored according to density and phase. The color hue is 
chosen from the standard color wheel based on the phase angle — i log(tp k / \ip x |). The corresponding 
eigenvalue approximation is Ax = 6.469449. The solution exhibits several vortices arranged in a 
regular pattern, as expected for a rotating BEC. 
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Figure 5.2. Plot of norm of residual and absolute error as function of iteration count (a) 
and ODE time t (b). The horizontal axis is linear in the left panel, and logarithmic in the right 
panel. In the final iterations, the time step is h max = 10 4 . The vector /j, denotes the residual, i.e., 
fk = P(v k )v k - A(y k )v k . 
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Figure 5.3. The approximations \ipk{ x >y)\ f or different iterates. The random starting value 
first turns into the general shape of the solution in ~ 10 iterations, and in the remaining iterations, 
only the vortices are modified. The contour levels were selected as lCT 10 , 10~ 8 , 10" 6 , 1CT 4 , 10" 3 , 
10 -2,5 , 10~ 2 , 10 -1 ' 75 , 10 -1 - 5 and 10 -1 ' 25 . The outer contours have the smaller values. 
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Figure 5.4. Estimation of the convergence factor from a calculation with h max = 250 and 
comparison with the theoretical convergence factor 7 (horizontal lines). The asymptotic region 
starts at iteration k RJ 95. The deviations for large k is due to numerical noise in the error estimate 
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Figure 5.5. Plot of o^ — A*, which quickly becomes small and negative. 



IEEE Transactions on Neural Networks, 1994. 

Appendix A. Derivation of local error and step-length for a variant 
of Rosenbrock-Euler. The error of the Rosenbrock-Euler method has been stud- 
ied in the context of Runge-Kutta methods (e.g. [HI Chapter IV. 7]). We need a 
more specialized result and will for completeness provide an error estimate for the 
Rosenbrock-Euler method in our setting. Consider the autonomous ODE 

y'(t) = f(y(t)) 

with ||y(0)|| — 1. Suppose the ODE has a structure such that the norm is an invariant, 
i.e., ||y(t)|| = 1 for all t > 0. Let y\ be one step of the Rosenbrock-Euler method. 
By using Taylor expansion, it is straightforward to show that the local error of the 
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Rosenbrock-Euler step is given by 



y 1 -y(h) = h 2 q+^-y'"(r) 





where 



(A.i) ?:=-^(/-/i/'(yo)))- 1 (/ + / i / / (yo)))/ / (yo)/(yo). 

and y'"(r) = (y'"(Ti)i,...,2/'"(r„)„) T , with ri,...,r„ e [0,h]. We approximate ||q|| w 
|||/'(yo)/(yo)|| an d neglect the final term which leads us to the following choice of 
step-length for a given error tolerance 



(A.2) h = 



2s 



\\r(m)f(mW 
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